Code
pacman::p_load(sp, sf, spNetwork, tmap, classInt, viridis, tidyverse, spatstat, sfdep, raster, maptools, arrow, lubridate)February 25, 2024
Human mobility, the movement of human beings in space and time, reflects the spatial-temporal characteristics of human behavior. With the advancement Information and Communication Technologies (ICT) especially smart phone, a large volume of data related to human mobility have been collected. By using appropriate GIS analysis methods, these data are potentially useful in supporting smart city planning and management.
In Singapore, one of the important source of data related to human mobility is from Land Transport Authority (LTA) DataMall. Two data sets related to human mobility are provided by the portal, they are: Passenger Volume by Origin Destination Train Stations and Passenger Volume by Origin Destination Bus Stops. One of the limitation of these data sets is that their location are biased to either bus stops or MRT/LRT stations. In 2020, another very interesting human mobility data set called Grab Posisi was released by GRAB, one of the largest shared taxi operator in South-east Asia. There are two data sets been released and one of them is for Singapore.
Geospatial analytics hold tremendous potential to address complex problems facing society. In this study, you are tasked to apply appropriate spatial point patterns analysis methods to discover the geographical and spatio-temporal distribution of Grab hailing services locations in Singapore.
The specific tasks of this take-home exercise are as follows:
Using appropriate function of sf and tidyverse, preparing the following geospatial data layer in sf tibble data.frames:
Grab taxi location points either by origins or destinations.
Road layer within Singapore excluding outer islands.
Singapore boundary layer excluding outer islands
Using the extracted data, derive traditional Kernel Density Estimation layers.
Using the extracted data, derive either Network Kernel Density Estimation (NKDE) or Temporal Network Kernel Density Estimation (TNKDE)
Using appropriate tmap functions, display the kernel density layers on openstreetmap of Singapore.
Describe the spatial patterns revealed by the kernel density maps.
For the purpose of this assignment, Grab-Posisi of Singapore will be used.
Road data set from OpenStreetMap of Geofabrik download server. The Malaysia, Singapore, and Brunei coverage should be downloaded.
Master Plan 2019 Subzone Boundary (No Sea) from Data.gov.sg.
The R packages used in this project are:
sf, a relatively new R package specially designed to import, manage and process vector-based geospatial data in R.
tidyverse for performing data science tasks such as importing, wrangling and visualising data.
spatstat, which has a wide range of useful functions for point pattern analysis. In this hands-on exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.
raster which reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this hands-on exercise, it will be used to convert image output generate by spatstat into raster format.
maptools which provides a set of tools for manipulating geographic data. In this hands-on exercise, we mainly use it to convert Spatial objects into ppp format of spatstat.
tmap which provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.
spNetwork, which provides functions to perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network. It also can be used to build spatial matrices (‘listw’ objects like in ‘spdep’ package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances.
rgdal, which provides bindings to the ‘Geospatial’ Data Abstraction Library (GDAL) (>= 1.11.4) and access to projection/transformation operations from the PROJ library. In this exercise, rgdal will be used to import geospatial data in R and store as sp objects.
sp, which provides classes and methods for dealing with spatial data in R. In this exercise, it will be used to manage SpatialPointsDataFrame and SpatiaLinesDataFrame, and for performing projection transformation.
Reading layer `MPSZ-2019' from data source
`C:\j00b00\IS415-GAA\Take-home_EX\Take-home_Ex01\data\Geospatial\MPSZ-2019'
using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84
Reading layer `gis_osm_roads_free_1' from data source
`C:\j00b00\IS415-GAA\Take-home_EX\Take-home_Ex01\data\Geospatial\malaysia-singapore-brunei-latest-free'
using driver `ESRI Shapefile'
Simple feature collection with 1759836 features and 10 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 99.66041 ymin: 0.8021131 xmax: 119.2601 ymax: 7.514393
Geodetic CRS: WGS 84
Before we proceed with our analysis, we need to preprocess the data to the correct format and clean it.
Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]
Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]
spatstat requires the analytical data in ppp object form. There is no direct way to convert a Spatial* classes into ppp object. We need to convert the Spatial classes* into Spatial object first.
class : SpatialPoints
features : 28000
extent : 3638.685, 50024.92, 25350.05, 49469.41 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
convert the spatial data into spatstat’s ppp object format.
Planar point pattern: 28000 points
Average intensity 2.502667e-05 points per square unit
Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units
Window: rectangle = [3638.69, 50024.92] x [25350.05, 49469.41] units
(46390 x 24120 units)
Window area = 1118810000 square units
check the duplication in a ppp object
When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.
Window: polygonal boundary
37 separate polygons (29 holes)
vertices area relative.area
polygon 1 71 5.63061e+03 8.47e-06
polygon 2 10 1.99717e+02 3.01e-07
polygon 3 12667 6.63014e+08 9.98e-01
polygon 4 (hole) 3 -3.41897e-05 -5.14e-14
polygon 5 (hole) 23 -1.99656e+01 -3.00e-08
polygon 6 (hole) 35 -1.38385e+02 -2.08e-07
polygon 7 (hole) 19 -4.39650e+00 -6.62e-09
polygon 8 (hole) 270 -1.21455e+03 -1.83e-06
polygon 9 (hole) 3 -4.95057e-02 -7.45e-11
polygon 10 (hole) 3 -3.65499e-03 -5.50e-12
polygon 11 (hole) 38 -7.79904e+03 -1.17e-05
polygon 12 (hole) 3 -5.99535e-04 -9.02e-13
polygon 13 (hole) 3 -3.04561e-04 -4.58e-13
polygon 14 (hole) 3 -7.43616e-06 -1.12e-14
polygon 15 (hole) 6 -8.37554e-01 -1.26e-09
polygon 16 (hole) 4 -2.86396e-01 -4.31e-10
polygon 17 (hole) 3 -1.81439e-04 -2.73e-13
polygon 18 (hole) 3 -8.68789e-04 -1.31e-12
polygon 19 (hole) 3 -4.46076e-04 -6.71e-13
polygon 20 (hole) 3 -3.39794e-04 -5.11e-13
polygon 21 (hole) 317 -5.11280e+04 -7.69e-05
polygon 22 (hole) 5 -2.92235e-04 -4.40e-13
polygon 23 (hole) 3 -4.52043e-05 -6.80e-14
polygon 24 (hole) 3 -3.90173e-05 -5.87e-14
polygon 25 (hole) 5 -2.44411e-04 -3.68e-13
polygon 26 (hole) 4 -2.18616e-04 -3.29e-13
polygon 27 (hole) 4 -4.28453e-01 -6.45e-10
polygon 28 (hole) 4 -2.54488e-04 -3.83e-13
polygon 29 (hole) 3 -9.59850e-05 -1.44e-13
polygon 30 (hole) 41 -4.01660e+04 -6.04e-05
polygon 31 (hole) 3 -4.14099e-04 -6.23e-13
polygon 32 (hole) 5 -3.64686e-02 -5.49e-11
polygon 33 30 2.80002e+04 4.21e-05
polygon 34 27 1.50315e+04 2.26e-05
polygon 35 285 1.61128e+06 2.42e-03
polygon 36 91 1.49663e+04 2.25e-05
polygon 37 71 8.18750e+03 1.23e-05
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units
(53270 x 28810 units)
Window area = 664597000 square units
Fraction of frame area: 0.433
Extract Grab origin point that are in main island in Singapore
Planar point pattern: 27821 points
Average intensity 4.186147e-05 points per square unit
Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units
Window: polygonal boundary
37 separate polygons (29 holes)
vertices area relative.area
polygon 1 71 5.63061e+03 8.47e-06
polygon 2 10 1.99717e+02 3.01e-07
polygon 3 12667 6.63014e+08 9.98e-01
polygon 4 (hole) 3 -3.41897e-05 -5.14e-14
polygon 5 (hole) 23 -1.99656e+01 -3.00e-08
polygon 6 (hole) 35 -1.38385e+02 -2.08e-07
polygon 7 (hole) 19 -4.39650e+00 -6.62e-09
polygon 8 (hole) 270 -1.21455e+03 -1.83e-06
polygon 9 (hole) 3 -4.95057e-02 -7.45e-11
polygon 10 (hole) 3 -3.65499e-03 -5.50e-12
polygon 11 (hole) 38 -7.79904e+03 -1.17e-05
polygon 12 (hole) 3 -5.99535e-04 -9.02e-13
polygon 13 (hole) 3 -3.04561e-04 -4.58e-13
polygon 14 (hole) 3 -7.43616e-06 -1.12e-14
polygon 15 (hole) 6 -8.37554e-01 -1.26e-09
polygon 16 (hole) 4 -2.86396e-01 -4.31e-10
polygon 17 (hole) 3 -1.81439e-04 -2.73e-13
polygon 18 (hole) 3 -8.68789e-04 -1.31e-12
polygon 19 (hole) 3 -4.46076e-04 -6.71e-13
polygon 20 (hole) 3 -3.39794e-04 -5.11e-13
polygon 21 (hole) 317 -5.11280e+04 -7.69e-05
polygon 22 (hole) 5 -2.92235e-04 -4.40e-13
polygon 23 (hole) 3 -4.52043e-05 -6.80e-14
polygon 24 (hole) 3 -3.90173e-05 -5.87e-14
polygon 25 (hole) 5 -2.44411e-04 -3.68e-13
polygon 26 (hole) 4 -2.18616e-04 -3.29e-13
polygon 27 (hole) 4 -4.28453e-01 -6.45e-10
polygon 28 (hole) 4 -2.54488e-04 -3.83e-13
polygon 29 (hole) 3 -9.59850e-05 -1.44e-13
polygon 30 (hole) 41 -4.01660e+04 -6.04e-05
polygon 31 (hole) 3 -4.14099e-04 -6.23e-13
polygon 32 (hole) 5 -3.64686e-02 -5.49e-11
polygon 33 30 2.80002e+04 4.21e-05
polygon 34 27 1.50315e+04 2.26e-05
polygon 35 285 1.61128e+06 2.42e-03
polygon 36 91 1.49663e+04 2.25e-05
polygon 37 71 8.18750e+03 1.23e-05
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units
(53270 x 28810 units)
Window area = 664597000 square units
Fraction of frame area: 0.433
Deriving kernel density estimation (KDE) layer for visualising and exploring the intensity of point processes
Re-run density() using the resale data set and plot the output kde map.

Properties of kde_originSG_bw_raster RasterLayer
class : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.4162063, 0.2250614 (x, y)
extent : 2.667538, 55.94194, 21.44847, 50.25633 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
source : memory
names : v
values : -3.709752e-13, 1953.626 (min, max)
Plot these four study areas and the locations of the origin points




kde_origin_pg.bw = density(origin_pg_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
kde_origin_tm.bw = density(origin_tm_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
kde_origin_ck.bw = density(origin_ck_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
kde_origin_jw.bw = density(origin_jw_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")







Drop redundant columns
Tampines
Choa Chu Kang
Jurong West
pg_samples$density <- pg_samples$density*1000
tm_samples$density <- tm_samples$density*1000
ck_samples$density <- ck_samples$density*1000
jw_samples$density <- jw_samples$density*1000
pg_lixels$density <- pg_lixels$density*1000
tm_lixels$density <- tm_lixels$density*1000
ck_lixels$density <- ck_lixels$density*1000
jw_lixels$density <- jw_lixels$density*1000Punggol
Observation : According to the Punggol map, it appears that Punggol Way has the highest density. This might be due to the limited convenience of public transportation in the area, leading residents to frequently utilize Grab services for transportation.
Tampines
Observation: According to the Tampines map, Changi Airport stands out as one of the areas with the highest density. This is logical, considering that individuals arriving in Singapore often rely on Grab services, especially since they may have luggage with them.
Choa Chu Kang
Observation : The Choa Chu Kang indicates that the central area od Choa Chu Kang, particularly in proximity to the MRT station, exhibits the highest density.
Jurong West
Observation : As per the Jurong West Map, the area with the highest density is located along Jurong West Avenue 2.